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Abstract. We present numerical investigations of the short-time dynamics at criticality in the ID Potts 
model with power-law decaying interactions of the form l/r^^" . The scaling properties of the magneti- 
zation, autocorrelation function and time correlations of the magnetization are studied. The dynamical 
critical exponents 9' and z are derived in the cases q = 2 and g = 3 for several values of the parameter a 
belonging to the nontrivial critical regime. 
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1 Introduction 

Short-time dynamics (STD) in systems quenched to criti- 
cality has attracted considerable attention in the last decade 
due to the appealing fact that systems even in the early pe- 
riod of relaxation to equilibrium exhibit universal scaling 
properties which involve both static and dynamic critical 
exponents [112] . The interest in this phenomenon exists at 
different levels. From a practical point of view, it offers 
a useful numerical tool for calculating both dynamic and 
static critical properties where the critical slowing down 
is turned into advantage. From a fundamental point of 
view, it opened a series of questions of current interest 
from the universal amplitudes to the universality of the 
fluctuation-dissipation ratio [|3j in a wider context of age- 
ing phenomena in pure systems [4J . One of the first points 
of conceptual interest was the emergence of a new indepen- 
dent universal dynamical exponent describing the initial 
increase of the magnetization in this regime [T , but related 
also to the persistence probability of the global order pa- 
rameter [S] . Since the STD was formulated in the context 
of the dynamical renormalization group (RG) and the new 
exponent evaluated within the e-expansion [IJ it has been 
further investigated, mostly numerically, in a variety of 
models in two and three dimensions for equilibrium phase 
transitions |6|7|8|9j and also for out-of-cquilibrium ones 

m- 

Quite a few studies were carried out on models with 
long-range (LR) interactions. The RG approach of Janssen 
et al. [1] was extended to the case of power-law decaying 
interactions of the form r~'^~'^ in the same continuous n- 
vector model [11], in the random Ising model [12j, and 



in the kinetic spherical model |13|14j . Studies of STD at 
criticality in discrete models with LR interactions, where 
such an approach does not apply, are still absent. Numer- 
ical "advantage" is there rather reduced due to the fast 
relaxation in the presence of LR interactions. 

In this paper we present the first and preliminary nu- 
merical study of the ID LR Potts model, useful as a para- 
digm that comprises different universality classes obtained 
by variation of the number of states q. We show that, in 
spite of the difficulties of the numerical approach in the 
LR case, the scaling properties characteristic for the STD 
may be well reproduced with a reasonable numerical ef- 
fort and derive the two dynamical critical exponents in 
the wide extent of the range-parameter a for two different 
universality classes. 

The outline of the paper is as follows. In Section [2] we 
give an overview of the model and basic STD properties 
considered in the paper, followed by the details of our 
numerical approach. The Section [3]contains the results for 
two special cases of the Potts model: q = 2, corresponding 
to the Ising model, which is compared to the previous RG 
results, and q — 3, where the new results are derived in 
the regime where the transition is of the second order. The 
conclusion is given in Section 2) 



2 Model and short-time dynamics approach 

We consider the ID Potts model defined by the Hamilto- 
nian 

^^ = -E ^4tT7'5.... ' (1) 



i<3 
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where J > 0, denotes a g-state Potts spin at the site 
i, 5 is the Kronecker symbol and the summation is over 
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all the pairs of the system. Hereafter J = fcs = 1 is 
used. As is well known [151, for < a < I the model 
([1]) has a phase transition at nonzero temperature for all 
q. Only a few exact results are available for its equilibrium 
critical behavior, but the model was studied in detail by 
several approximate methods [16117118] . It has a rather 
complicated phase diagram in the {q, a) plane, involving 
similar variety of critical regimes, that is encountered in 
the {q, D) plane of the same model with short-range (SR) 
interactions. This gives the additional motivation to ex- 
amine also the dynamical scaling properties in the STD 
regime depending on q and a. 

In the present work we are interested in two special 
cases, q = 2 and q = 3 in the range of parameter a 
corresponding to the nontrivial (non mean- field (MF)) 
critical regime, where the initial slip of the magnetiza- 
tion can be observed. For q — 2 this is accomplished for 
0.5 < CT < 1 [Ej. In the latter case, q — S, which belongs 
to a different universality class, this region is restrained to 
adq = 3) < (7 < 1, where <Tc(q) > 0.5 denotes the point 
of the onset of the first-order phase transition, occurring 
for q > 2 and known only approximately |20l21j . For these 
two cases we shall study the nonequilibrium evolution to 
criticality in early times of several quantities, magneti- 
zation, autocorrelation function and time correlations of 
the magnetization. Let us first briefly remind their scaling 
properties in the STD regime and explain their implemen- 
tation to the model ([T]). 



2.1 STD approach 

As shown by Janssen et al [T], if the system is brought 
out of equilibrium by a quench from high temperature to 
criticality, and left to evolve following the nonconservative 
dynamics of Model A (in the sense of reference |22]), then, 
during the early stage of relaxation it will display univer- 
sal scaling properties characterized by the static exponents 
and the new universal dynamic exponent. Consequently, 
in the system of size L after a quench from high temper- 
ature to the critical region in the presence of small initial 
magnetization mo, the magnetization will obey the scaling 
relation 

M{t, T, L, mo) = b-^^/''M{t/b',b^/''T,L/b, 6^«mo), (2) 

where r = {T—Tc)/Tc, 6 is a scaling factor and /3, ly are the 
static critical exponents. Besides the dynamical exponent 
z, the scaling involves a new exponent xq as the anomalous 
dimension of the initial magnetization mo. 

At criticality (r = 0), and for i ^, equation ^ may 
be reduced to 



M{t,mo) = t 



Xq/z 



ma) 



-z/xo 




(3) 

but larger 



For early times satisfying t tx ~ 
than the microscopic time tmicro, the r. h. s. can be ex- 
panded giving the power-law increase of the magnetization 
known as the initial slip. 



M (t) - moi^ , 



(4) 



with 6' = xq/z — [3 1 {vz). The magnetization in the model 
^ is defined in a standard way 



M{t) = (Mi(t)) 



{q-l)L 



E 



(5) 

where 1 denotes the preferential direction among q possi- 
ble Potts states a. The brackets (...) denote the average 
over initial conditions and random force. 

During the short time after the quench, the correlation 
length is small compared to the system size, and the expo- 
nent 9' can be derived directly from the power law ^ by 
performing simulations on the chain of a single large size 
and averaging over a great number of independent runs. 

In the absence of the initial magnetization (mo — 0), 
equation Q gives the scaling relation for the fc-th moment 
of the magnetization. 



M'^''\t,L) = b-''^'''M'^''\tlb\L/b). 



(6) 



In early times, when ^{t) <^ L, the second moment also 
displays a power-law behavior. 



(7) 



which can be used to derive the anomalous dimension of 
the order parameter P/v^ or the dynamical exponent z 
directly from the single large chain. To this purpose we 
use the alternative definition of the order parameter 



M,(t) 



{q-l)L 



(8) 



and the moments of magnetization are obtained as the 
average 

M«(t) = (A/,^(i)). (9) 

Equation ([8]) describes the absolute value of the magne- 
tization and allows us to apply the scaling relation ([6]) 
already to the first momentum, and obtain 



{Mx{t)) - t'^d/2-f)/u)/z^ 



(10) 



The autocorrelation function of the local order param- 
eter is defined in a standard way and also obeys the power- 
law form 



A{t) 



q 



{q-l)L 



E 



-\/z 



(11) 

depending on the combination of both dynamical expo- 
nents \/ z — d/ z — 9' . 

For the calculation of the exponent 9' we shall use 
another quantity which represents the autocorrelation of 
the global order parameter. It was shown by Tome and de 
Oliveira lOj that the time correlation of magnetization 
defined as 

Q{t) = {Mi{0)Mi{t)) (12) 
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also exhibits the initial increase of the power-law form 

Q{t) ~ t"' (13) 

even in absence of the imposed initial magnetization. (No- 
tice that in Equation the definition ^ of magneti- 
zation should be used and not its absolute value.) For 
numerical calculations of the exponent 9\ Equation (|13p 
has a technical advantage compared to the expression (g]), 
where the runs should be performed first for several val- 
ues of the initial magnetization mo and than the extrap- 
olation to the limit mg — > taken in order to obtain the 
exponent 9' . In return, however, the fluctuations are more 
pronounced for Q{t) and its calculation requires better 
statistics. 



2.2 Numerical calculations 



Table 1. Number of steps elapsed before the correlation length 
(|14p has reached half of its maximum value for q — 2. 



a\L 


100 


400 


1000 


3 000 


0.9 


2 


13 


40 


168 


0.8 


1 


7 


20 


54 


0.7 


1 


4 


10 


30 


0.6 





3 


6 


15 



were performed with chains of 3 000 sites. All the quan- 
tities were averaged over 200 000 to 350 000 independent 
runs. Larger numbers of independent runs were used for 
smaller values of a, where the fluctuations are more pro- 
nounced. Finally, in the ESS approach small sizes ranging 
from L = 100 to L = 400 were compared. 



Monte Carlo simulations were done on flnite chains with 
periodic boundary conditions by using simple Metropo- 
lis dynamics. The system was quenched from a random 
configuration (high-temperature state) to criticality. 

Unlike the earlier studies for the short-range Potts 
model in 2D, where the critical temperatures are known 
exactly, in the LR case only the approximate results are 
available. Satisfactory results for series of different values 
of a were obtained by the finite-range scaling (FRS) ap- 
proach [23], cluster mean-field approach [24], or Monte 
Carlo calculations [U]. In the present study we use the 
values for Tc obtained by the FRS [i6j . 

Two approaches were examined - a direct derivation 
of exponents and a derivation from the finite-size scal- 
ing (ESS). In the former approach where the exponents 
are calculated using a single large system, the correlation 
length has to be small compared to the system size dur- 
ing times which are taken into account in the evaluation 
of the power laws. Due to LR interactions, the correla- 
tion length increases much faster than in systems with 
short-range interactions. For illustration we supply here 
a rough estimate of the increase of the correlation length 
^(t) calculated from the second moment of the spin-spin 
correlation function [25j at the instant t, 

where the correlation fmiction C{l,t) is given by 

^(m) = (^(e(w..(*)-^))- (15) 

Summation in equation p4p runs only up to L/2 because 
of the periodic boundary conditions. As shown in Table [H 
the correlation length increases very rapidly indeed, espe- 
cially for lower values of a. Values for q = 3 are similar. 
Consequently, in order to reach sufficiently long time in- 
tervals in the power-law regime, all the direct calculations 



3 Results 

Systematic calculations in cases q — 2 and q = 3 were 
performed for four characteristic values of parameter of 
range a = 0.6, 0.7, 0.8, 0.9. 

Increasing of a by moving away from the ME regime 
up to the limits of relevance of long-range interactions has 
similar effect in this ID model as leaving the ME regime 
by lowering dimensionality down to the lower critical di- 
mensionality in its SR analogue, and we expect to observe 
similar features. One of them is dependence of dynamical 
exponents 9' and z on a. 

The above choice of a allows to cover evenly the non- 
trivial critical regime 0.5 < cr < 1 for q = 2. In the case 
q = 3 it covers both first- and second-order transition 
regimes, but the detailed analysis is focused on the region 
where the second-order phase transition is expected. 



3.1 Case q = 2 

3.1.1 Time correlations of the magnetization and the 
exponent 9' 

The dynamical exponent 9' , which in the SR analogue in- 
creases with decreasing of dimensionality [1] , in the present 
LR case should increase with cr, which is also in agreement 
with the RG results [TT] . 

The principal quantity that we used to derive the ex- 
ponent 9' is the function Q{t) p2| . A summary graph of 
our numerical simulations for the selected values of a is 
presented in Figure [1] One observes that the microscopic 
time tjnicro is short and the linear behavior on the log- 
arithmic scale is established immediately after the first 
2-3 steps. The linear regime in the log-log scale becomes 
shorter as a decreases, but the size L = 3 000 is sufficiently 
long for an accurate evaluation of the exponent 9' by a fit 
to equation (|13p . which deteriorates only for the lowest a 
considered. 
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Table 3. The exponent 6' calculated from the magnetization 
for q = 2, a = 0.9 and L = 3 000 for several values of initial 
magnetization mo with the linear extrapolation to mo = 0. 



mo 


0.1 


0.05 


0.01 


mo 


e'imo) 


.187 ± .005 


.196 ± .006 


.202 ± .008 


.204 ± .009 



using this approach. Just for illustration, we present the 
cr = 0.9 data in Table[3l limiting ourselves to a very rough 
estimation. 

The result is consistent with the one cited in Table 
[21 Improving the accuracy would imply performing the 
calculations on several smaller initial values mo, each of 
them requiring the same amount of numerical effort spent 
for the calculation of Q{t). 



Fig. 1. The time correlations of the magnetization for q = 2 
and the selected values of cr (L = 3 000). 



3.1.2 Magnetization and exponent z 



The errors in present results are not easy to estimate, 
because they may be introduced by several sources: insuf- 
ficient statistics, arbitrariness in the selection of the lin- 
ear segment of the plot, or using the approximate values 
for Tc- The error bars given in tables cover the first two 
sources and could be systematically reduced by increasing 
the number of independent runs and the size of the chains 
considered. Yet, we stress that the third one cannot be 
estimated directly. 

The obtained values for 6' are presented in Table [2] 
compared to the RG results by the two-loop e expansion 
of reference |11| . Since the accuracy of our results improves 



Table 2. Dynamical exponents for q — 2 compared to the RG 
results of reference |Tl] rounded up to 4 digits. 



a 


e' 






z 




ZRG 


Xjz 


0.9 


0.212 ± 


.005 


0.3346 


1.18 ± 


.04 


0.9532 


0.635 ± .004 


0.8 


0.188 ± 


.004 


0.2587 


0.96 ± 


.04 


0.8340 


0.85 ±.01 


0.7 


0.137 ± 


.006 


0.1733 


0.81 ± 


.01 


0.7174 


1.136 ± .02 


0.6 


0.07 ± 


.01 


0.0821 


0.70 ± 


.01 


0.6052 


1.47 ± .02 



with increasing cr, one may conclude that the RG results 
are overestimated due to the insufficiency of the two-loop 
expansion in that regime. Similar overestimation of 6' was 
observed in the SR case, where the same RG e-expansion 
in the SR limit [11] gives e.g. & = 0.131 and 9' = 0.356 for 
d — 2> and d — 2 respectively, while the MC simulations 
give respectively 9' = 0.104 [55] and 9' = 0.191 [HT]. 

A more standard way to calculate the exponent 9' is 
from the initial slip of the magnetization given by equation 
In the present problem we find it less advantageous 
both for precision and for the numerical effort needed. For 
this reason we do not proceed with the systematic analysis 



As discussed earlier in Section 2.2, the dynamical expo- 
nent z is expected to increase with cr, since the relaxation 
becomes slower with decreasing range of interactions. 

The exponent z was calculated from the magnetiza- 
tion using equation pi7)) . The log-log plot of the simulation 
data is illustrated in Figure [D for the case a — 0.8. The 



Int 



Fig. 2. Magnetization data for q = 2, a — O.l 
the linear extrapolation (dash-dotted line). 



(diamonds) and 



values of z presented in Table [2] were obtained by substi- 
tuting into equation (llOp the exact value for the anoma- 
lous dimension of the order parameter P/v, which is equal 
to (1 — cr)/2 [23. As cr increases our results become signif- 
icantly larger than those obtained by the e expansion |llj . 
Again, we may attribute this discrepancy to an underes- 
timation of the RG results by the two loop expansion and 
may observe similar behavior in the SR case, where, for 
the Ising model, the (4 — d)-expansion to the second order 
[28] gives z = 2.013 and z = 2.052 for d = 3 and d = 2 
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respectively, while the best MC calculations give z close 
to 2.04 g for d=3, and z==2.1667 ^ for d = 2. 

3.1.3 Autocorrelation function 

The example of simulations of the autocorrelation function 
for a — 0.8 is illustrated by the log-log plots in Figure 
[31 The power-law fit to equation (llip gives the exponent 



FigureElgives the scaling fit for the case q — 2 and a — 0.9. 
By Cimaa: wc deuotc the saturation value that attains 
according to the expression p4)) . It is proportional to the 
size L in the limit of large L. As in the previous example, 
the scaling is performed with the same value of z as given 
in Table H 

The agreement in both cases is very good. Neverthe- 
less, these fits are generally less accurate than direct cal- 
culations from systems of large sizes. 



< 

B .4.5 



Int 



Fig. 3. MC data for autocorrelation function for q = 2, a = 
0.8, L = 3 000 (circles) and the linear extrapolation (dotted 
line) . 

A/z presented in Table [21 Although not suitable for the 
calculation of the exponent 9' when the exponent z is not 
known with sufficient precision, the values for A/z were 
used for a check of independent calculations of & and z. 
Within given error bars, the agreement is obtained. 



3.1.4 Finite-size scaling 

An alternative way of evaluating the exponent z is to per- 
form the simulations on several small systems of different 
sizes and apply FSS by using the overlapping fits jJOJ . To 
this purpose one may consider the magnetization, Binder's 
fourth-order cumulant, but also the correlation length de- 
fined by equation (fT4)) . 

We illustrate two such fits, involving sizes L = 100, 200 
and 400, for the magnetization and the correlation length 
in Figures [H and [5l 

The fit for the magnetization was performed by apply- 
ing the scaling relation ^ to the magnetization defined 
by ([8]). An example for q — 2 and <t = 0.8 is given in 
Figure [H The magnetization is rescaled by using the ex- 
act value for (3 /v. The time axis is rescaled by using the 
earlier calculated value of the exponent z cited in Table [21 

The scaling fit may also be applied directly to the cor- 
relation length defined by equation , since at the crit- 
icality it should scale as 



3.2 Case q — 3 

For the three-state Potts model we expect to obtain dif- 
ferent dynamical exponents. Also, we should be able to 
distinguish, depending of a, two regimes, corresponding 
to the first- and second-order phase transition. 

The calculations were performed along the same lines 
and with similar parameters as for the preceding case, 
since the increase of the correlation length with time is 
very similar to that for q — 2. 

Yet, the microscopic time period tmicro was found to 
be larger by several steps than the one for q = 2. This 
property is clearly seen in Figure [6l In spite of such be- 




Ut)=L f{t/L'). 



(16) 



Int 

Fig. 6. Comparison of tmicro for the cases q — 2 and q — 3 
on the example of the function Q{t) (equation (|12|) ') for a — 
0.9, L — 3 000. The straight lines are the linear extrapolations 
performed in the scaling regime. The plot for q — 3 was shifted 
by 0.5 in the y direction in order to display the two plots on 
the same graph. 



havior, owing to the fact that the analysis for q = 3 is 
limited to larger values of a (as explained later), it was 
sufficient to use the same size L — 3 000 as in the q — 2 
case, but the statistics had to be increased systematically 
up to the 350 000 independent runs. 

In Figure [71 we present the results for the time corre- 
lation function of the magnetization, Q(t), for the same 
values of a as in the previous case q — 2. The behavior is 
qualitatively different from the one in Figure [H An initial 
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= 400 



Fig. 4. Scaled magnetization M^{t, L) for L = 100, 200, 400 (g = 2, cr = 0.8). 




Fig. 5. Scaled correlation length for L = 100, 200,400 (g = 2, cr = 0.9). 
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Fig. 7. The time correlation of the magnetization for g = 3 
and the selected values of a (L = 3 000). 



this connection, it is important to notice, that STD has 
proven as an efficient approach in cases involving short- 
range interactions for studying both first-order phase tran- 
sitions }32|33l34 ] and a the tricritical point [35]. A more 
detailed study of these issues in the present LR model re- 
quires a separate study [36j . We mention here only, that 
our preliminary results for Q{t) on a finer scale of a lo- 
cate the change of regime between 0.72 and 0.74, which 
is in agreement with most recent estimates |21| that give 
Gc{q = 3) = 0.72(1). Here we shall limit the scope to the 
second-order transition regime analyzing further only the 
behavior for a = 0.8 and a = 0.9. 

The dynamical exponent 6' derived from the log-log 
plot of the function Q{t) is given in Table H) It strongly 

Table 4. The dynamical exponents for g = 3. 

a 9' z ZFSS ^/z 

0.9 0.120 ±.004 1.21 ±.01 1.26 ±.04 0.704 ±.008 
0.8 0.058 ± .004 1.01 ± .004 1.02 ± .04 0.935 ± .006 



increase is observed only for a = 0.8 and a = 0.9, while 
the change of behavior for lower values of a announces the 
expected onset of the first-order transition regime, which 
lies between a = 0.7 and a = 0.8. The onset of the first- 
order phase transition in the present model is a challeng- 
ing question in itself, since the position of the tricritical 
point ac{q) is still not known with precision [18|31|21j . In 



decreases as the first-order regime approaches. Compared 
to the Ising case, the exponents 6' for q = 3 turn out to be 
significantly lower. This is similar to what was observed 
for the 2D SR Potts model, where 0' = 0.191 for q = 2, 
while 0' = 0.075 for g = 3 [7]. 



Katarina Uzelac et al.: Short-time dynamics in the ID long-range Potts model 



7 



In Table [5] we also present the alternative derivation 
of the exponent 9' by investigating the initial slip of the 
magnetization for the case a = 0.9. A rough linear ex- 



Table 5. The exponent 9' derived from the magnetization in 
the case g = 3, a = 0.9 (L = 3 000) for several values of the 
initial magnetization mo with the extrapolation to mo = 0. 



mo 


0.1 


0.05 


0.02 


mo 


6l'(mo) 


.133 ± .003 


.118 ± .004 


.112 ± .006 


.106 ± .09 



trapolation to toq = gives a slightly smaller value for 9' 
than the one cited in the Table ID In the same time the 
precision of the calculations from the magnetization was 
considerably lower. 

Within this approach, one can observe one feature com- 
mon to earlier numerical calculations for the 2D short- 
range Potts model [7], that 9' {mo) converge to the limit 
mo from different sides for q = 2 and q = 3, which 
was there attributed to the opposite positions of the re- 
lated fixed points. 

The results for the exponent z obtained from equation 
(fTO|) using the same procedure as in the case q — 2 are 
presented in the third column of Table [4] and illustrated 
in Figure [8] for a = 0.8. The values of z are slightly larger 
than for the Ising case, similar as it was obtained for the 
2D SR Potts model ^7^. 

For want of prior results for dynamical exponents in 
the case g = 3, we also applied the FSS to the magne- 
tization ^ and the correlation length (fT4|) by perform- 
ing independent evaluations of the exponent z by using 
the collapsing fits. The results are included in Table [4] for 
comparison. 

As in previous case, the results for X/ z obtained from 
the autocorrelation function (jlip (cited in the last column 
of Table [3]) agree with the independently calculated values 
of 0' and z within the accuracy limits. 



4 Conclusion 

We presented a numerical study of scaling properties re- 
lated to the short-time dynamics at criticality in the ID 
LR Potts model. Based on the analysis of several physical 
quantities, we showed that in spite of the fast relaxation in 
presence of the LR interactions, the STD scaling regime 
can be observed numerically and dynamic critical expo- 
nents evaluated with satisfactory accuracy. We focused 
here on studying the problem in larger range of a, but the 
accuracy of each individual result may still be improved 
with reasonable numerical effort. 

The dynamical exponents 9' and z were evaluated in 
the cases q = 2 and g = 3, for several values of a belonging 
to the nontrivial critical regime. The exponents are found 
to differ for the two cases and depend on a, in similar way 



T 




J , \ 1 \ , I , 

2 2.5 3 3.5 

Int 

Fig. 8. MC data for the magnetization (|10p in the case g = 3, 
a = 0.8, L = 3 000 (diamonds) and the linear extrapolation 
(line). 

they depend on dimensionality in the SR analogue of this 
model. 

For the Ising case, the comparison could be made with 
the existing RG results. A fair agreement for values of tr 
close to the MF border (e = 2cr — 1 ^ 1) is obtained, 
but the discrepancy reaches far beyond the estimated er- 
ror bars as a increases, which could be attributed to the 
shortcomings of the e-expansion. Our results are in favor 
of significantly smaller increase of Q' and larger increase 
of z with decreasing range of interactions. 

For g = 3, new values for the exponents & and z were 
obtained in, more restrained, second-order phase transi- 
tion regime. The value of the exponent z is found to be 
slightly larger than the one for g = 2, while increasing 
number of Potts states had larger impact on the critical 
exponent 9' which is appreciably smaller and tends to van- 
ish as the first-order transition regime approaches. We also 
found the change in the behavior of the time correlations 
of the magnetization as the first-order transition sets in 
with lowering of the parameter g. 

Besides the onset of the first-order transition regime in 
this model which is already a subject of a separate study, 
a number of issues remain to be examined further, such as 
the possible effects of different dynamics, or a complemen- 
tary analysis of the exponent describing the persistence 
probability of the global order-parameter at criticality. 

This work was supported by the Croatian Ministry of Science, 
Education and Sports through grant No. 035-0000000-3187. 
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